function [Q,grad,se,p1mat] = func_mixed_logit_current(theta,y,X,draws,sumidx,...
    XC,drawsC,sumidxC,seoption)

beta = theta(1:size(X,2));
sigma = theta(size(X,2)+1:end);

sumdraws = zeros(size(draws,1),size(draws,3));
sumdrawsC = zeros(size(drawsC,1),size(drawsC,3));
for k=1:size(draws,2)
   sumdraws = sumdraws + sigma(k)*squeeze(draws(:,k,:));
   sumdrawsC = sumdrawsC + sigma(k)*squeeze(drawsC(:,k,:));
end

umat = (X*beta)*ones(1,size(draws,3)) + sumdraws;
expumatC = exp((XC*beta)*ones(1,size(drawsC,3)) + sumdrawsC);


umatCfull = zeros(max(sumidx),size(sumdrawsC,2));
umatCfull(sumidxC,:) = expumatC;
umatCfull = umatCfull(sumidx,:);

XCfull = zeros(max(sumidx),size(XC,2));
XCfull(sumidxC,:) = XC;

XCfull = XCfull(sumidx,:);

drawsCfull = zeros(max(sumidx),size(drawsC,2),size(drawsC,3));
drawsCfull(sumidxC,:,:) = drawsC;

drawsCfull = drawsCfull(sumidx,:,:);


p1mat = max(min(exp(umat)./(1+umatCfull+exp(umat)),.9999999),.00000001);
p0mat = 1-p1mat;

pCmat = umatCfull./(1+umatCfull+exp(umat));

ymat = y*ones(1,size(draws,3));

probmat = p1mat.*ymat+p0mat.*(1-ymat);

Lmat = zeros(max(sumidx),size(draws,3));

for s=1:size(draws,3)
    Lmat(:,s) = accumarray(sumidx,probmat(:,s),[],@prod);
end

countidx = accumarray(sumidx,ones(size(sumidx)),[],@sum);

Lmatmean = mean(Lmat,2);

Q = -mean(log(Lmatmean(countidx>0)));

inner = (2*ymat-1).*(1-probmat).*Lmat(sumidx,:);

innerC = ((p1mat-ymat)./p0mat).*pCmat.*Lmat(sumidx,:);

grad = zeros(size(theta));

for b=1:length(beta)
    tempgrad = accumarray(sumidx,X(:,b).*mean(inner,2),[],@sum)./Lmatmean;
    tempgrad = tempgrad + accumarray(sumidx,XCfull(:,b).*mean(innerC,2),[],@sum)./Lmatmean;
    grad(b) = -mean(tempgrad(countidx>0));
end

for b=length(beta)+1:length(theta)
    tempgrad = accumarray(sumidx,mean(squeeze(draws(:,b-length(beta),:)).*inner,2),[],@sum)./Lmatmean;
    tempgrad = tempgrad+accumarray(sumidx,mean(squeeze(drawsCfull(:,b-length(beta),:)).*innerC,2),[],@sum)./Lmatmean;
    grad(b) = -mean(tempgrad(countidx>0));
end

if seoption==0
se = [];
else
varcov = zeros(size(theta,1));  

fderiv = zeros(max(sumidx),size(theta,1));

for b=1:length(beta)
    fderiv(:,b) = -accumarray(sumidx,X(:,b).*mean(inner,2),[],@sum)./Lmatmean;
    fderiv(:,b) = fderiv(:,b)-accumarray(sumidx,XCfull(:,b).*mean(innerC,2),[],@sum)./Lmatmean;
end

for b=length(beta)+1:length(theta)
    fderiv(:,b) = -accumarray(sumidx,mean(squeeze(draws(:,b-length(beta),:)).*inner,2),[],@sum)./Lmatmean;
    fderiv(:,b) = fderiv(:,b)-accumarray(sumidx,mean(squeeze(drawsCfull(:,b-length(beta),:)).*innerC,2),[],@sum)./Lmatmean;
end

for b=1:size(theta,1)
    for c=1:size(theta,1)
        varcov(b,c) = mean(fderiv(countidx>0,b).*fderiv(countidx>0,c));
    end
end

varcov = inv(varcov)/sum(countidx>0);

se = sqrt(diag(varcov));
end


